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I. INTRODUCTION 

Thanks to recent advancess in radioactive ion beam 
technology, we are now in the process of exploring the 
very limits of nuclear binding, namely those regions of 
the periodic chart in the neighborhood of the particle drip 
lines Several new structure features have already 

been uncovered in these studies, including the neutron 
halo, and others have been predicted. 

In contrast to stable nuclei within or near the valley of 
beta stability, a proper theoretical description of weakly- 
bound systems requires a very careful treatment of the 
asymptotic part of the nucleonic density. This is partic- 
ularly true in the description of pairing correlations near 
the neutron drip line, for which the correct asymptotic 
properties of quasiparticle wave functions and of one- 
particle and pairing densities is essential. In the frame- 
work of the mean-field approach, the best way to achieve 
such a description is to use the Hartree-Fock-Bogoliubov 
(HFB) theory in coordinate-space representation [^-|7|. 

Such an approach presents serious difficulties, how- 
ever, when applied to deformed nuclei. On the one hand, 
for finite-range interactions the technical and numerical 
problems arising when a two-dimensional mesh of spatial 
points is used are so involved that reliable self-consistent 
calculations in coordinate space should not be expected 
soon. On the other hand, for zero-range interactions ex- 
isting approaches |^,^ are able to include only a fairly 
limited pairing phase space. The main complication in 
solving the HFB equations in coordinate space is that the 
HFB spectrum is unbounded from below, so that meth- 
ods based on a variational search for eigenstates cannot 
be easily implemented. Because of this and other difh- 
culties, one has to look for alternative solutions. 

In principle, such an alternative solution is well known 
in the form of the configurational representation. In this 
approach, the system of partial differential HFB equa- 
tions are solved by expanding the nucleon quasiparticle 



wave functions in an appropriate complete set of single- 
particle wave functions. In many applications, an expan- 
sion of the HFB wave function in a large harmonic oscil- 
lator (HO) basis of spherical or axial symmetry provides 
a satisfactory level of accuracy. For nuclei at the drip 
lines, however, expansion in an oscillator basis converges 
much too slowly to describe the physics of continuum 
states 1^], which play a critical role in the description 
of weakly-bound systems. Oscillator expansions produce 
wave functions that decrease too steeply in the asymp- 
totic region at large distances from the center of the nu- 
cleus. As a result, the calculated densities, especially in 
the pairing channel, are too small in the outer region 
and do not reflect correctly the pairing correlations of 
such nuclei. 

In two recent works 
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a new transformed har- 
monic oscillator (THO) basis, based on a unitary trans- 
formation of the spherical HO basis, was discussed. This 
new basis derives from the standard oscillator basis by 
a local-scaling point coordinate transformation [p^-p^, 
with the precise form dictated by the desired asymptotic 
behavior of the densities. The transformation preserves 
many useful properties of the HO wave functions. Using 
the new basis, characteristics of weakly-bound orbitals 
for a square- well potential were analyzed and the ground- 
state properties of some spherical nuclei were calculated 
in the framework of the energy den sity functional ap- 
proach |ll). It was demonstrated in ||l^ that configura- 
tional calculations using the THO basis present a promis- 
ing alternative to algorithms that are being developed for 
coordinate-space solution of the HFB equations. 

In the present work, we develop the THO basis for use 
in HFB equations of axially-deformed weakly-bound nu- 
clei. Our main goal here is to present and test these new 
theoretical methods. As specific applications, we repeat 
previous calculations performed for the chain of Mg iso- 
topes H), but for different effective interactions, and then 
report a preliminary study of light, neutron- rich nuclei 



1 



near the drip line. Extensive calculations throughout the 
mass table, together with a more detailed analysis of the 
pairing interaction, will be presented in a future publica- 
tion. 

The structure of the paper is the following. The THO 
basis for deformed nuclei is introduced in Sec. O. In 



Sec. [II we present an outline of the HFB theory and 
discuss several features of particular relevance to our in- 



vestigation. Results of calculations are given in Sec. IV 
and conclusions are presented in Sec. 0. 



II. TRANSFORMED HARMONIC OSCILLATOR 
BASIS 

In this section, we introduce a generalized class of 
local-scaling point transformations, which in principle 
act differently in the three Cartesian directions. Next, 
we apply this transformation to the three-dimensional 
Cartesian HO wave functions and derive the correspond- 
ing properties of the local densities. 



A. Local-scaling point transformations 

Suppose {ipa(r)} represents a complete set of or- 
thonormal single-particle wave functions depending on 
the spatial coordinate r. (To simplify the presentation, 
we suppress the spin and isospin labels here.) Then, one 
can introduce a local-scaling point transformation (LST) 
of the three dimensional vector space, which is a gen- 
eralization of the analogous spherically-symmetric LST 
namely 



X — > x' = x'{x,y,z) ^ ffr,{r), 
y — > y' = V'{x, y, z) = f /y(r), 
z > z' = z'{x,y,z) = ^f,Xr), 



where r—^Jx"^ +2/^4- 2^. 

The LST functions fk{r), k=x, y, or z_ 



mathematical properties ensuring that (2.1 



(2.1) 



should have 
is a valid 



invertible transformation of the three-dimensional space. 
In particular, /fc(r) should be monotonic functions of r 
such that 



/fe(0) = and /fe(oo) 



(2.2) 



and should lead to a non- vanishing Jacobian of the LST 
(P), i.e.. 



D : 



d{x\y',z') 
d{x,y, z) 

fyj z I y jxjyj 



fxfyfz y fxfyfz + Z fxfyfz 



(2.3) 



where primes denote derivatives with respect to r. 

When we apply the LST ( 2.1) to the set of wave func- 
tions ipair), we obtain another set of single-particle wave 
functions 



y, z) = D^'\^{^J,{r), ^fy{r), f /,(r)). (2.4) 



Due to the factor D^/^ entering Eq.(|2j), the LST of wave 
functions is unitary and the new wave functions tpaif) are 
automatically orthonormal, i.e. , {ipa\ipi3) — {'Pa\^i3)—Sai3- 
Summarizing, the LST (2.1) generates from a given 
complete set of orthonormal single-particle wave func- 
tions another orthonor mal and complete set of single- 
particle wave functions (|]J) depending on three almost- 
arbitrary scalar LST functions fk{r). The freedom in the 
choice of /fe(r) provides great flexibility in the new set 
{?Act(r)}, and this opens up the possibility of improving 
on undesired properties of the initial set. This is the mo- 
tivation for the present study in which we use the LST to 
modify the incorrect asymptotic properties of deformed 
HO wave functions. 



B. Transformed harmonic oscillator wave functions 

The anisotropic three-dimensional HO potential with 
three different oscillator lengths 



has the form 
U{r) 



Lk = T- 



(x" 



mujk 



2m\Li Ll Li 



(2.5) 



(2.6) 



Its eigenstates, the separable HO single-particle wave 
functions 

'fair) = (pnAx)fny{y)'PnAz), (2.7) 

have a Gaussian asymptotic behavior at large distances. 



fair oo) ~ exp 



2 V -t'x -^y -^2 



(2.8) 



Applying the LST (2.1) to these wave functions leads 
to the so-called THO single-particle wave functions (|]J), 

^„(r) = D'/' {fUr)) (^„„ (f /,(r)) (f /.(r)), 

(2.9) 

whose asymptotic behavior is 



i'air oo) ^ exp 



1 h'f'x . y'fy . 



2 \ £2^2 



-I- 



(2.10) 



This suggests that we choose the LST functions to satisfy 
the asymptotic conditions 



fk{r) = 



r for small r, 

LkV^KT for large r. 



(2.11) 
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With such a choice, the THO wave functions at small r 
are identical to the HO wave functions (note that with 



(1 



11) one obtains D=l at small r), while at large r they 
have the correct exponential and spherical asymptotic 
behavior, 



(2.12) 



C. Parametrization of the LST functions 

In principle, we could use the flexibility of having three 
different LST functions fk{r) and three different oscilla- 
tor lengths Lk of the original deformed HO basis to tailor 
the LST transformation to the shape of the deformed nu- 
cleus under investigation. However, for large HO bases 
(in the present study we include HO states up to 20 major 
shells), the dependence of the total energy on the basis 
deformation is very weak, so that minimization of the 
total energy with respect to the three oscillator lengths 
Lk is ill-conditioned (see discussion and examples given 
in Ref. |jl^). Therefore, in this study, we use a spherical 
HO basis depending on a single common oscillator length 



Lx — Ly — Lz = Lq — — 
bo 



(2.13) 



With such a choice, it is natural to set the three LST 
fimctions /fe(r) equal to one another. 



Ur) = !y{r)=fz{r)^f{r). 



(2.14) 



This allows us to use exactly the same LST function 
/(r) a s in the previous studies ||l^,|ll| . Under conditions 
( 2.14 ), the Jacobian (2^) assumes the simpler form 



D . 



dix',y\z') ^ .r{r)p{r) 
d{x,y,z) r2 



(2.15) 



The parametrization of the LST function /(r) used in 
Refs. [PUnl was of the form 



(2.16) 



with the dimensionless universal function F of the di- 
mensionless variable TZ defined as 



'7^(l-^a7^2 



,1/3 



for TZ < c . 



F(7^): 



d-2 



n 



-do + diTZ+dLlnTZ for 7^ > c 



(2.17) 



Two different formulae can be obtained for the function 
F{TZ), one for 7^ < c and one for TZ > c. Imposing 
the condition that the function should be continuous at 



the matching radius c and that it should have continuous 
first, second, third, and fourth derivatives leads to the 
following requirements for the constants d_2, d^i, do, 
di, and dL- 



rf_2=3^c4(243 
d_i=-8^c3(81 
do 



40507 - 
- 12427 



di 
dL 



89107^ + 76027^ -|- 22757^^), 
' 27457^ + 234O73 -I- 7007^^), 
=2y^c2(12157 + 27907^-^241573 + 72874 -I- (486 
+64807 + I43IO72 + I2I8O73 + 36407-*) In c), 
= |ylc(243 + 243O7 + 52II72 + 438O73 + 13007^), 
=-4y^c2(243 + 324O7 + 715572 + 6O9O73 + 18207^), 

(2.18) 



where 7=00^ and A "'^=81(1 -f 7) 
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In this way, the 



LST function /(r) is guaranteed to be very smooth, while 
still depending on only three parameters, Lq, a and c. 

From (2.17), we see that asymptotically the function 
F{TZ- ^oo)^ ^/ diTZ. Thus, the LST function obeys condi- 
tions ( |2.11 ) provided that the parameters satisfy 

K=^. (2.19) 

Two different approaches can be used in calculations. 
One possibility is to minimize the total energy with re- 
spect to Lq, a, and c, obtaining as output the energet- 
ically optimal value of the decay constant k. Alterna- 
tively, for a given choice of k, we could eliminate one of 
the three parameters and minimize the total energy with 
respect to the other two. The actu al pro cedure used in 



our calculations is described in Sec. IV A 



D. Axially deformed harmonic oscillator 

In the present study, we restrict our HFB analysis to 
shapes having axial symmetry. For this purpose, we use 
HO wave functions in cylindrical coordinates, z, p, and 



X = pcosip , 
y = psinip , 



(2.20) 



which allows us to separate the HFB equations into 
blocks with good projection Q of the angular momentum 
on the symmetry axis. [Note that the use of cylindrical 
coordinates is in dependent of working with equal oscilla- 
tor lengths ( 2.13| ) .] Since the use of a cylindrical HO basis 
is by now a standard technique (see, e.g., Ref. [0), we 
give here only the information pertaining to constructing 
the cylindrical THO states. 

The cylindrical HO basis wave functions are given ex- 
plicitly by 



iPa{z,p,ip,S,t) = (p„Jz).^™'(p)— =Xme('S)XmtW 



2tt 



(2.21) 
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where the spin s and isospin t degrees of freedom are 
shown explicitly, and Up are the number of nodes along 
z and p directions, respectively, while mi and iris are the 
components of the orbital angular momentum and the 
spin along the symmetry axis. The only conserved quan- 
tum numbers in this case are the total angular momen- 
tum projection il=mi+ms and the parity 7r=(— )"^+™' 



In the axially deformed case, the general LST {2A) acts 
only on the cylindrical coordinates z and p and takes the 
form 



P — ' P' ^ P'{p.z) = f 

z — > z' = z'{p,z) = Uzir), 



with the corresponding Jacobian given by 
D 



_ dix',y',z') _ P^f'pfpf. + z^flr. 



d{x,y,z) 

Finally, the axial THO wave functions are 



(2.23) 



(2.24) 



The assumption of a single oscillator length (see Sect. 



II C ) that we make in our calculations translates in the 
axial case to 



Lp — Lz 



n 



/p(r) = /.(r)^/(r). 



(2.25) 



(2.26) 



and the Jacobian (2.23) reduces to expression (p.l5|) 



E. THO and Gauss integration formulae 



(1 



At first glance, the THO wave functions (2.9) and 
.24 ) loo k m uch m ore co mplicated than their HO coun- 
terparts ( ^.TD and (2.21). In particular, in contrast to 
the HO wave functions, the THO wave functions are not 
separable either in the x, y, and z Cartesian coordinates 
or in the p and z axial coordinates. Due to the presence 
of the Jacobian factor and the r-dependence of the LST 
functions, the local-scaling transformation mixes the x, y 
and z coordinates and the p and z coordinates. Neverthe- 
less, as we now proceed to show, the THO wave functions 
are readily tractable in any configurational self-consistent 
calculation. Indeed, the modifications required to trans- 
form a code from the HO to the THO basis are minor. 

One of the properties of the HO basis that makes it 
so useful is the high accuracy that can be achieved when 
calculating matrix elements using Gauss-Hermite and/or 
Gauss-Laguerre integration formulae |17|. This feature 
has been exploited frequently in various mean-field nu- 
clear structure calculations (see, e.g., Refs. p^JT^,|l5t). 



To illustrate how the same methods can be applied in 
the THO basis, we focus on the specific example of a 
diagonal matrix element of a spin and isospin indepen- 
dent potential function V. This matrix element can be 
expressed in the axial HO representation as 



oo oo 



ifclVl^^) = / dz pdpV{z,p)^l{z)^^f{p), 



-oo 



(2.27) 



(2.22) and in the THO representation as 



(V'al^l^a) 



dz p dpV{z, p) 



X D{z,p) i^fzir)) (£/^(^)). (2.28) 



The way to calculate the second matrix element ( 2.2S ) 
is by first transforming to the p' and z' variables (2.22). 
This absorbs the Jacobian D{z,p) and leads to an inte- 
gral o ver HO wave functions that is almost identical to 
(p.27D, namely 



OO OO 



(V'al^^lV'o) = J dz' J p' dp'V{z{z',p'),p{z',p')) 
-oo 

X ^ljz')^^'\p'). (2.29) 



The only complication in numerically carrying out 
the integral ( ^.29| ) involves determining the inverse LST 
transformations z=z{z', p') and p=p{z' , p') to be inserted 
into the known function V{z, p). But this only has to be 
done once, and, moreover, if Gauss quadratures are used 
to evaluate the integrals, the inverse transformation only 
has to be known at a finite number of Gauss-quadrature 
nodes. 

Generalization of the above approach to include differ- 
ential operators, as will often arise in THO basis configu- 
rational calculations, is fairly straightforward. Such inte- 
grals can be done by first transforming derivatives d/dz 
and d/dp into derivatives d/dz' and d/dp', and then per- 
forming the integrations in the variables z' and p' over 
ordinary HO wave functions (see the next section). 



F. THO and local densities 

In calculations using the Skyrme force, or in any other 
calculation that relies on the local density approximation, 
we can simplify the THO methodology of Sec. II E even 
further. Indeed, suppose the mean-field calculation in 
question relies on knowing the density matrix paa' in the 
THO basis. Then the spatial nonlocal density can be 
expressed as 



P{ri,r2) 



>a (r^Paa'C' (^2) , 



(2.30) 
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and the corresponding standard local densities as 



p{r) = p{r,r) 
T(r) 

jkir) 



k—x^y,z 

1 

2i 



where 



(2.31a) 
(2.31b) 

(2.31c) 



(2.32) 



for 1=1 or 2, an d k=x, y, or z. To simplify the nota- 
tion in Eq. ( ^.30| ), we have neglected the spin and isospin 
degrees of freedom and, conse quent ly, have shown only 
the spin- independent densities ( 2.31 ). Analogous formu- 
lae for the spin-dependent densities s, T, and Jm (l^ are 
straightforward. 

A direct calculation of the derivatives in E qs 
[after inserting the THO wave functions (2.9) or ( 



2.31) 



2.24) 



into the nonlocal density matrix (2.30)] is prohibitively 
difficult. Fortunately, nothing of the sort is necessary. 
It is enough to note that the densities ( 2.31| ) serve 
almost uniquely to define the central, spin-orbit, and 
effective-mass terms of the mean-field Hamiltonian (see, 
e.g., Refs. ]T9| , p^ ), and that these terms are in turn used 
t o cal culate matrix elements thro ugh in tegrals of the type 

have to be effec- 



(2.29). Therefore, the densities (2.31 



y', z' (the Gauss- 



tively known only at selected points x' , 
quadrature nodes) of the inverse LST. 

Towards this end, we insert the THO wave functions 
into the nonlocal density ( 2.30| ), which gives 



p{r,,r2) = D'/\r,)D'/\r,)p'{r[ir,),r'^{r2)), (2.33) 



with 



(2.34) 



The density matrix p'{r'i,r2) is a standard object ex- 
pressed in terms of ordinary HO wave functions, and it 
can be calculated using methods that are employed in 
any code that works in the HO basis. Likewise, the cor- 
responding local densities 



p'(r')-p'(r',r') 



Vi^^'v(^)'p'(r-;,r',) 



(2.35a) 
(2.35b) 



can be calculated without any reference to the THO ba- 
sis. The only complication is that now we have to cal- 
culate the complete kinetic energ y tenso r density 
( 2.35b| ), while finally only its trace (2.31b) is needed. In- 
serting expression (2.33) into (2.31), and expressing the 
differential operators~(2.32) as 



V 



(i) 



for 



^mk q 



dVk 

we obtain that 

p{r{r'))=Dp'{r'), 

T{r{r')) = DY^D'^D'^rUr' 



(2.36) 

(2.37) 
(2.38a) 



iE[Vfc^]^r[v;,p'(r-') 



--D-'[VDfp{r'), 



Jk 



Mr'))^DY,DT3'm{r'). 



(2.38b) 
(2.38c) 



To use formulae ( 2.3^ ), we must calculate the Jacobi 
matrix and its determinant D at points r{r')] how- 
ever, this need be done only once for all iterations. On 
the other hand, no inverse LST nee ds to be performed for 
the densities, because expressions ( 2.38 ) give directly the 
values of the local densities at the inverse LST po ints, as 
required in matrix- element integrals of the type ( ^.29 ). 



III. HARTREE-FOCK-BOGOLIUBOV THEORY 

Hartree-Fock-Bogoliubov (HFB) theory is based 
on the Ritz variational principle applied to the many- 
fermion Hamiltonian, 

-ff = E + E ^aa'i3i3'al^al,a,3'ai3, (3.1) 

aa' aa'f3f3' 

with trial functions in the form of a quasiparticle vacuum. 
The resulting HFB equations can be written in matrix 
form as 

where En are the quasiparticle energies, A is the chemical 
potential, and the matrices h = t + T and A are defined 
by the matrix elements of the two-body interaction 



013' 



(3.3) 



Pfi'p and Kpfji being the density matrix and pairing ten- 
sor, respectively. HFB theory is by now a standard tool in 
nuclear structure calculations, and we refer the reader to 
Ref. [|o| for details. Below we discuss several features of 
the formalism that are especially pertinent to the present 
application, namely canonical states, the pairing phase 
space, and those quantities that dictate the stability of a 
nucleus with respect to two-neutron emission. 
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A. Canonical states 

Canonical states are defined as the states that diag- 
onaU ze th e HFB one-body density matrix p{ri,r2) of 
Eq. (^!30|) , i.e., 



(3.4) 



where, due to the Pauli principle, the canonical occupa- 
tion numbers vf obey the condition < vf < 1. 

For self-consistent solutions, the canonical occupation 
numbers vf are determined by the diagonal matrix ele- 
ments hii and A^^ of the particle-hole (p-h) and particle- 
particle (p-p) Hamiltonians in the canonical basis via the 
following BCS-like equation po| : 



A 



2E,, 



where 



A2 



(3.5) 



(3.6) 



The chemical potential A is determined from the particle 
number condition 

i n 

where iV„ denote the norms of the lower HFB wave func- 



tions of Eq. (B^), i.e.. 



(3.8) 



In the canonical representation, the average (proton or 
neutron) pairing gap A |^ is given by the average value 
of Aii in the corresponding (proton or neutron) canonical 
states. 



(3.9) 



wher e N is the number of nucleons of that type (see 
Eq. (P)). 

Whenever infinite complete single-particle bases are 
used in configurational calculations, one may freely 
expa nd the upper and lower HFB wave functions of 
Eq. (3^) (the quasiparticle wave functions), as well as 
the standard eigenstates of the p-h Hamiltonian h, in the 
canonical basis. These expansions are often extremely 
slowly converging, however, and any truncation of the 
basis typically induces large errors. Therefore, in prac- 
tice, when working with finite bases, one should not ex- 
pand quasiparticle wave functions and the single-particle 
eigenstates of h in the canonical basis. The reason is very 
simple, and it stems from different asymptotic properties 
of these objects. As discussed in Ref. ||], the quasiparti- 
cle spectrum and wave functions are partly discrete and 



localized and partly continuous and asymptotically os- 
cillating, respectively. These properties are completely 
analogous to properties of the eigenstates of h, which 
are also discrete and localized (for negative eigenener- 
gies) or continuous and oscillating (for positive eigenen- 
ergies). On the other hand, the proper ties of eigenvalues 
and eigenstates of the density matrix (3.4) are very dif- 
ferent, namely the entire spectrum is discrete and all the 
wave functions are localized. Therefore, even if formally 
the set of canonical states is complete, it is extremely 
difficult to expand any oscillating wave function in this 
basis. 

These considerations make it clear that the optimum 
way of solving the HFB equations is by using the coor- 
dinate representation, in which the various asymptotic 
properties are in a natural way correctly fulfilled. This 
technique is widely used when spherical symmetry is 
imposed; then one only has to solve systems of one- 
dimensional differential equations, which is an easy task. 
On the other hand, the case of axial symmetry requires 
solving two-dimensional equations, and that of triaxial 
shapes requires working with a three-dimensional prob- 
lem. None of these two latter cases has up to now been 
effectively solved in coordinate space, although work on 
the axial solutions is in progress pl|] . 

Therefore, without having access to coordinate- 
representation solutions, we are obliged to use methods 
based on a configurational expansion. In this respect, 
one may clearly distinguish two classes of finite single- 
particle bases, each of which aims at a reasonable solu- 
tion of the HFB equations (3.2). One uses a truncated 
basis composed of eigenstates of ft. [22 This basis is 



partly composed of discrete localized states and partly of 
discretized continuum and oscillating states. Technically 
it is very difficult to include many continuum states in 
the basis, especially when triaxial deformations are al- 
lowed. In practice, Refs. included states up to 
several MeV into the continuum. Such a small phase 
space is certainly insufficient to describe spatial proper- 
ties of nuclear densities at large distances, although some 
ground-state properties, like total binding energies, will 
be at most weakly affected. 

The second uses a truncated infinite discrete basis. 
The most common of course is the HO basis, which 
has been used in numerous HFB calculations, espe- 
cially those employing the Gogny effective interaction 
(see, e.g., Refs. p^-|2^), and in Hartree-Bogoliubov 
calculations based on a relativistic Lagrangian (see, 
e.g., Refs. p^ , p8t ). Because it uses a basis with a similar 
structure to the canonical basis (infinite and discrete), 
this approach can be viewed as aiming at the best pos- 
sible approximation to the canonical states and not the 
quasiparticle states. In this sense, the amplitudes Un and 
Vn that appear in Eq. ( |3.2| ) should be considered more 
as expansion coefficients of quasiparticle states in a basis 
similar to the canonical basis than as quasiparticle wave 
functions themselves. 

Our approach, which we discuss in greater detail be- 
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low, belongs to the second class. The THO basis defined 
and described in Sec. || is a model that aims at an op- 
timal description of the canonical states. Therefore, in 
the following we adapt properties of the THO b asis, and 
in particular the value of the decay constant k (2.12), to 
the asymptotic properties of canonical states. In fact, the 
unique decay constant of all THO basis states is exactly 
the desired property of canonical states. As discussed 
in Ref . , the asymptotic properties of the most impor- 
tant canonical states (those having average energies close 
to the Fermi energy) are governed by a common unique 
decay constant. 



'2m(i;,nin- A) 



(3.10) 



where -Bmin is the lowest quasiparticle energy En- This 
should be contrasted with decay constants associated 
with the eigenstates of which are all different and de- 
pend on the single-particle eigenenergies. 



B. The cut-off procedure 

HFB calculations in configurational representation in- 
variably require a truncation of the single-particle basis 
and a truncation in the number of quasiparticle states. 
The latter is usually realized by defining a cut-off quasi- 
particle energy Smax and then including quasiparticle 
states only up to this value. When the finite-range Gogny 
force is used both in the p-p and p-h channels, the cut-off 
energy -Emax has numerical significance only. In contrast, 
HFB calculations based on Skyrme forces in the p-h and 
p-p channels, as well as any other calculations based on 
a zero-range force in the p-p channel 



y*(r,r') ==yo'5(r-r') 



(3.11) 



require a finite space of states. This is because, for 
any value of the coupling constant Vb, they give diver- 
gent energies with increasing -Emax (see the discussion in 
Ref. 0). 

The choice of an appropriate cut-off procedure has 
been discussed in the case of coordinate-space HFB cal- 
culations for spherical nuclei 1^. It was demonstrated 
there that one must sum up contributions from all states 
close in quasiparticle energy to the bound particle states 
to obtain correct density matrices in the HFB method. 
Since the bound particle states are associated with quasi- 
particle energies smaller than the absolute value D of the 
depth of the effective potential well, one had to take the 
cut-off energy i?max comparable to D. 

In the case of deformed HFB calculations, and espe- 
cially when performing configurational HFB calculations, 
it is difficult to look for the depth of the effective poten- 
tial well in each subspace. Thus, an alternative cri- 
terion with respect to the above cut-off procedure used 
in spherical calculations is needed. For this purpose, we 



have adopted the following procedure (see Appendix B of 
[^). After each iteration, performed with a given chemi- 
cal potential A, we calculate an auxiliary spectrum e„ and 
pairing gaps A„ by using for each quasiparticle state the 
BCS-like formulae, 



En 

or equivalently 



(e-„-A)^ + A2, 

1 _ en - A 

2 2 En 



(1 - 2Nn)En 



An - 2En^Nn{l~Nn). 



(3.12a) 
(3.12b) 



(3.13a) 
(3.13b) 



Then, in the next iteration, we readjust the proton and 
neutron chemical potentials to obtain the correct val- 
ues of the proton and neutron particle numbers ( [3.7| ), 
whe re agai n Nn is calculated for the equivalent spectrum, 
Eq. (3.12b). Due to the similarity between the equivalent 
spectrum e„ and the single-particle energies, we are tak- 
ing into account only those quasiparticle states for which 



fin — 



(3.14) 



where eniax>0 is a parameter defining the amount of 
the positive-energy phase space taken into account. At 
the same time, since all hole-like quasiparticle states, 
Nn<} /2, have negative values of e„ ( 3.13a ), condition 
(|3.14 ) quarantees that they are all taken into account. 
In this way, we have a global cut-off prescription inde- 
pendent of , which fulfills the requirement of taking 
into account the positive-energy phase space as well as 
all quasiparticle states up to the highest hole-like quasi- 
particle energy. 



C. Two- neutron separation energies and Fermi 
energies 

A particular thrust of our analysis will be to iden- 
tify the location of the two-neutron drip line. The 
self-consistent HFB variational procedure produces two 
quantities that provide information of relevance. One is 
the two-neutron separation energy, S2m defined as the 
difference between the HFB energy for the A*"— 2 and N 
neutron systems (with the same proton number) and the 
other is the Fermi energy, A„. 

The two-neutron separation energy provides "global" 
information on the total Q-value corresponding to a hy- 
pothetical simultaneous transfer of two neutrons into the 
A^— 2 ground state, leading to the ground state of the 
nucleus with N neutrons. The Q- value includes infor- 
mation on all differences in the ground-state properties 
of both nuclei, like pairing, deformation, configuration, 
etc. Whenever this Q-value becomes negative, the win- 
dow for the spontaneous and simultaneous emission of 
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two neutrons opens up, and the nucleus with N neutrons 
is formally beyond the two-neutron drip line. 

The Fermi energy, on the other hand, gives "local" in- 
formation on the stability of the given nucleus at a given 
pairing intensity, deformation, and configuration. Within 
the HFB theory, the sign of the Fermi energy dictates the 
localization properties of the HFB wave function; it is lo- 
calized if A„<0 and unlocalized (i.e., behaves asymptoti- 
cally as a plane wave) if A„>0. Thus, within the HFB ap- 
proach, nuclei with Xn>0 spontaneously emit neutrons, 
while those with A„<0 do not emit neutrons, irrespec- 
tive of the available Q- values for the real emission. As 
such, we must take into account all solutions with A„<0 
in discussing our self-consistent HFB results. 

We will indeed see examples in Sec. IV in which the 
nucleus has a negative two-neutron separation energy, so 
that it is formally beyond the two-neutron drip line, but 
nevertheless is localized and does not spontaneously spill 
off neutrons. 



IV. RESULTS 

In this section, we present the results of several sets of 
HFB calculations performed in the axial-deformed THO 
basis. All the calculations were carried out using the 
Skyrme interaction SLy4 |29j, which has recently been 
adjusted to the properties of stable nuclei, neutron-rich 
nuclei and neutron matter. This force has a proven record 
in deformed- mean-field calculations ||30|-|3^,p|,^3|-|36| , in- 
cluding calculations of rotational properties in nuclei 
P7H41|. At the same time, it reproduces the masses of 
spherical nuclei with an accuracy similar to several other 
Skyrme forces. 

Below we review our choice of the various parameters 
that define our calculations and present several tests of 
the THO approach. Then we present results obtained for 
the Mg isotopes and for light nuclei at the two-neutron 
drip line. A more extensive set of calculations will be 
presented in a future study, where we shall also explore 
in detail the influence of the type of pairing force on the 
properties of drip-line nuclei. 



A. Parameters and numerical details of the 
calculation 

In all of th e cal culations reported here, we use a contact 
interaction ( 3.11 ) in the p-p channel, which leads to vol- 
ume pairing correlations 0|. Following the discussion of 
Sec. IIIB , the pairing phase space has been defined by a 
cut-off energy (see Eq. (3.14)) of emax=30MeV. This con- 
stitutes a very safe limit, for which all convergence prop- 
erties are well satisfied (see the discussion in Refs. ||7,p2|). 
Wit hin th is phase space, the pairing strength Vq (see 
Eq. (3.11)) has been adjusted in a manner analogous to 



the prescription used in Re f. [|43| , namely so that the av- 
erage neutron pairing gap (|3.9D for ^^'^Sn equals the ex- 
perimental value of A„=1.245MeV. The resulting value 
is Vo=— 206MeVfm'^. As demonstrated in the Appendix 
of Ref. 0, changes in the cut-off parameter Cmax, lead- 
ing to a renormalization of the pairing strength Vq, can 
be safely disregarded when compared to all other uncer- 
tainties in the methods used to extrapolate to unknown 
nuclei. 

Although our axially-deformed HFB-I-THO code is 
able to work with arbitrary axial oscillator lengths Lp 
and Lz, we have used in these calculations a spherical 
basis defined b y a single common oscillator length Lq 
( p.25D (see Sec. |l^). When optimizing the THO basis 
parameters Lq, a, and c (to minimize the total energy), 
we invariably find that for weakly-bound nuclei the re- 
sulting exponential decay constant ( 2.19 ) is very close to 
that given by the HFB estimate (3.1C). Based on this 
observation, we have chosen to eliminate the THO pa- 
rameter a and to fix it in such a way that the basis decay 
constant (2.19), at the self-consistent solution, is equal 
to the HFB decay constant (3.10). In this way, we only 
have two variational parameters in our calculations, Lq 
and c. The minimizations were carried out independently 
for each nucleus. When describing each specific applica- 
tion, we will indicate the number of shells included both 
in the minimization that determines the LST parameters 
(^siD ^^'^ final calculations (A^sh)- 

All Gauss integrations were performed with 22 nodes 
in the p direction and 24 nodes in the z direction (due to 
the reflection symmetry assumed with respect to the x-y 
plane, only 12 nodes for z>0 were effectively needed). 



B. Tests of the method 

As the first test of the method, we considered doubly- 
magic nuclei. Such nuclei are known to be spheri- 
cal and thus amenable to reliable calculation using the 
coordinate-space HFB code Q. By studying the extent 
to which our code is able to reproduce the coordinate- 
space results (referred to in subsequent discussion as ex- 
act), we can assess the method. [We should note here 
that HFB in fact reduces to HE in doubly-magic nuclei, 
since all pairing correlations vanish] . 

These calculations were carried out both for nuclei 
along the beta-stability line and for the very neutron-rich 
nucleus ^*0. Some discussion of this latter nucleus is in 
order here. ^^O is known experimentally to be unbound 
Q,^, but is predicted to be bound in most mean- field 
calculations Due to rapid changes of the single- 

particle energies with neutron number, shell- model cal- 
culations are able to explain the sudden decrease of 
separation energies that occurs in the chain of oxygen 
isotopes and renders ^^O and ^^O unbound. This effect 
seems to require modifications to the effective interac- 
tions currently in use in mean-field studies of light nuclei. 
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Nevertheless, it is common to use ^^O as a testing ground 
of mean-field calculations near the neutron drip line, be- 
cause according to the standard magic-number sequence 
it is doubly magic and because it is located (in typi- 
cal mean-field calculations) just before the two-neutron 
drip line. This is the philosphy underlying our inclusion 
of ^^O. For comparison, the configurational calculations 
were performed both in the HO and THO bases. To as- 
sess the convergence of the results in the two cases, we 
varied the number of HO major shells included, consid- 
ering Nsh—8, 12, 16, and 20. For a given number of the 
major shells, we minimized the total HF energies with 
respect to the basis parameters, Lq for the HO basis, 
and Lq and c for the THO basis, so here N^-^'^—Nsh- We 
also tested our HO axial-basis results obtained at any 
given A'sh with those available from Cartesian-basis cal- 
culations ijl^ and the results agreed perfectly. Lastly, 
for the THO basis, wc compared with the calculations 
of Rcf . fill ] , where spherical symmetry was imposed, and 
obtained identical results. 

As expected, for nuclei within the valley of beta sta- 
bility the HO and THO results are close to one another 
and, furthermore, coincide with the exact HFB (HF) re- 
sults. The situation is quite different for the neutron-rich 
nucleus ^*0, for which the calculations indicate the pres- 
ence of a significant neutron skin. In Fig. ^, we present 
the HO and THO results for the total energy and for the 
proton and neutron rms radii as functions of A^sh- For 
each of the calculated observables, the exact results are 
shown as a straight line as a function of Nsh- Clearly, 
when we increase the number of major shells, both the 
HO and THO results for the total energy and for the 
proton rms radius converge to the exact HFB values. In 
contrast, the HO neutron rms radius still differs from the 
exact value, even at A'sh=20, while the THO basis gives 
the correct result. 

An explanation of this difference becomes clear when 
looking at Fig. ||, in which we compare (in logarithmic 
scale) the HO and THO neutron densities with those 
from the exact HFB calculations. The HO neutron den- 
sity fails to reproduce the correct asymptotic behavior 
at large distances (see also the discussion in Ref. 0). 
The THO density, on the other hand, shows perfect 
agreement with the exact HFB density. There is a dif- 
ference, of course, near and beyond the box boundary 
(^box=20 fm is used in the coordinate HFB calculations). 
The coordinate-space density rapidly falls to zero at the 
boundary, while the THO density continues with the cor- 
rect exponential shape out to infinite distances. 

It is clear that the rather small numerical discrepancy 
between the HO and THO neutron rms radii (Fig. |^) 
does not reflect the seriousness of the error in neutron 
densities that arises when using the HO basis. It is also 
obvious that observables which do not strongly depend on 
neutron densities at large distances, like the total energy 
or proton radii, are fairly well reproduced in standard 
HO calculations. On the other hand, observables that 
do depend on densities in the outer region, most notably 



pairing correlations Q], require the correct asymptotic 
behavior provided by the THO basis. 

Encouraged by the excellent results in spherical nuclei, 
where a comparison with reliable coordinate-space calcu- 
lations was possible, we next turned to deformed systems. 
Here, since no coordinate-space HFB results are avail- 
able, our tests were limited to a study of the convergence 
of results with increasing number of HO shells. [The ex- 
act results would be obtained in either the HO or THO 
expansion with a complete space, i.e., an infinite number 
of shells.] Whenever the number of HO shells used in 
the final HFB calculation was 12 or less, we determined 
the basis parameters with that same number of shells, 
N^^'^=Nsh- When the number of HO shells of the final 
calculations exceeded 12, however, we still determined 
the basis parameters with N^{^' = 12. 

In Fig. |, we show convergence results for the ground 
state of the weakly-bound deformed nucleus ^°Mg. The 
top three panels give the results for the total energies, the 
proton rms radii and the neutron rms radii, respectively. 
The fourth gives results for the /3 deformation, which is 
related to the quadrupole moment (Q) (Q — X^iLi '^^f ~ 



and the rms radius 



by 



(Q) 

5 A{r^) 



(4.1) 



The results obtained with A^sh=20 arc indicated in the 
figure by horizontal lines. Again, both bases yield very 
good convergence for the total energy and proton ra- 
dius. In contrast, noticeable differences between the HO 
and THO results can be seen for the deformation and 
neutron rms radius, and they persist to large values of 
Ash. Although these differences are small in magnitude, 
they are caused by a very large error in the HO neu- 
tron density distribution. This is illustrated in Fig. 
where we show the neutron densities calculated for the 
nearby ^^Mg nucleus. Every point in the figure corre- 
sponds to the value of the neutron density at a given 
Gauss-integration node. Since there are always several 
nodes near a sphere of the same radius r = 
there can be some scatter of points, corresponding to dif- 
ferent densities in different directions. This is especially 
true at small distances. At large distances, the scatter 
is greatly reduced and the densities exhibit to a good 
approximation spherical asymptotic behavior, exponen- 
tial in the case of the THO expansion and Gaussian in 
the case of the HO expansion. Note, however, that some 
scatter persists in the THO results out to large distances, 
suggesting that deformation effects are still present there. 
This is apparently reflecting the importance of deforma- 
tion of the least-bound orbitals. Clearly, the asymptotic 
properties of the HO and THO neutron densities are very 
different from one another, as they were in the spherical 
calculations (see Fig. 0). 
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C. Drip-line-to-drip- line calculations in Mg 

The chain of even-Z magnesium isotopes has been the 
subject of numerous recent theoretical analyses. The ex- 
treme interest in this isotopic chain is motivated by re- 
cent measurements in '^^Mg psj-pof, which show a larger- 
then-expected quadrupole coUectivity. Based on the rel- 
ativistic and non-relativistic mean-field approaches and 
on shell- model calculations (see Ref. ||5l[] for a review), it 
is now well documented that shape coexistence and con- 
figuration mixing occur in this N—20 nucleus. Moreover, 
recent advances in radioactive-ion-beam technology allow 
mass measurements of even heavier isotopes [^2| , giving 
hope that the neutron-drip line can be experimentally 
reached in the Z=12 chain 

In this section, we present results of an investiga- 
tion of the deformation properties of the even-even Mg 
isotopes from the proton-drip-line to the neutron-drip- 
line. Our results are complementary to those of recent 
Skyrmc+HFB calculations [||, in which the imaginary- 
time evolution method of finding eigenstates of the mean- 
field Hamiltonian h (see Sec. IV) was combined with a 
diagonalization of the HFB Hamiltonian within a rela- 
tively small set of these eigenstates. In that study, a com- 
plete set of results was given only for the SIII force and 
density-dependent pairing was used. Here, we present a 
complete set of results for the SLy4 force with a density- 
independent (volume) pairing interaction. These calcu- 
lations were carried out with Nsh—^O and N^^'=12 HO 
shells. 

In Fig. |, we plot the total HFB energies per nucleon 
E/A, the neutron chemical potentials A„, the neutron 
and proton deformation parameters, /?„ and Pp, the neu- 
tron, proton, and total quadrupole moments, Qm Qp, 
and Qt, the average neutron and proton pairing gaps, 
A„ and Ap, Eq. (|3.9|), and the pairing energies E^^^^. 
and -Epair for the magnesium isotopes as functions of the 
mass number A. Ground-state values are shown by full 
symbols connected by lines, while the isolated open sym- 
bols correspond to secondary minima of the deformation- 
energy curves. In the top panel of Fig. |^, we compare the 
results for the two-neutron separation energies S2n (open 
symbols) with those for the related quantity — 2A„ (full 
symbols), and in the bottom panel we show the neutron 
and proton rms radii. 

The lightest Mg isotope predicted by these calcula- 
tions to be bound against two-proton decay is ^°Mg. The 
heaviest bound against two- neutron decay, on the basis of 
having a positive two-neutron separation energy, is '^''Mg. 
On this basis, the position of the two- neutron drip line 
obtained within the HFB-|-SLy4 approach is identical to 
that obtained in the finite-range droplet model l54| , rela- 
tivistic mean field (RMF) ||5|], and HFB+SIII ||fcalcu- 
lations. The RMF approach with the NL-SH effective in- 
teraction |5^] predicts the two-neutron drip line at ^^Mg, 
and the relativistic Hartree-Bogoliubov (HB) approach 
with the NL3 effective interaction |M predicts it at or 



beyond ''^Mg. 

On the other hand, from Fig. ^, we see that both ''^Mg 
and ^''Mg, though having negative values of S2n, have 
(small) negative values of the Ferm i energy, A„ . Accord- 
ing to the discussion of Sec. IIIB, these nuclei, both of 
which exhibit oblate shapes, are bound against neutron 
emission. We will return to this point later. 

The most deformed nucleus of the isotope chain is ^"^Mg 
with almost the same neutron and proton deformations. 
At the other end of the chain, due to a large excess of neu- 
trons over protons, significant differences exist between 
the proton and neutron quadrupole moments. The on- 
set of large deformation in ■^^Mg causes a decrease of the 
neutron chemical potential A„ with respect to its value in 
■^'^Mg. This gives an additional binding of ^^Mg, and cor- 
respondly to an increase and decrease of the two-neutron 
separation energies S'2ri in "^^Mg and "^^Mg, respectively 
(see Fig. In experiment |^^, these changes are less 
pronounced and arrive two mass units earlier, giving rise 
to a small and large decrease of S2n in '^''Mg and "^^Mg, 
respectively. 

Concerning the ground-state deformation properties 
(full symbols connected by lines in Fig. the proton 
drip-line nucleus ■^"Mg displays a well defined spherical 
minimum {N—8 is a magic number). Then, there is 
a competition between prolate (^^'^''Mg and 36,38,40]y[g-j 
and oblate (^^''^"Mg) deformations, while ^^^-^^Mg are 
spherical. The last two localized isotopes (with nega- 
tive Fermi energies), '^^'''''Mg, display oblate deforma- 
tions. Secondary minima of the deformation energy 
curves (isolated symbols) exist for isotopes '^'^•'^'^•'^^Mg and 

36,38,40]yjg 

Non-zero proton pairing correlations are present at 
all spherical or oblate minima. However, at these 
shapes, tangible neutron pairing exist only in ^^'^''Mg and 
34,36,38]y|-g Moreover, for all nuclei with prolate ground- 
state shapes, i.e., in ^^'^'^Mg and 36,38,40Mg^ y^^^^-^ proton 
and neutron correlations are small or vanish altogether. 
These results are at variance with the Gogny-pairing HB 
calculations of Ref. , where non-zeropairing exists in 
all the heavy Mg isotopes. Also, in Ref. g], stronger pair- 
ing correlations were obtained for the zero-range density- 
dependent pairing force. However, in that study, the 
strength parameters were not adjusted to odd-even mass 
staggering but rather taken from high-spin calculations 
of superdeformed bands. Our results suggest that the 
pure HFB-pairing approach is not necessarily the best 
way to treat pairing correlations in the Mg isotopes, and 
approximate or exact particle-number projection should 
probably be employed. 

At this point, it is worth expanding a bit on the un- 
usual results for ''^'^^Mg. In these two isotopes, the 
solutions corresponding to prolate shapes are unsta- 
ble (A„>0), while those corresponding to oblate shapes 
continue to be bound, i.e., they have A„=— 0.253 and 
-0.092 MeV for A=A2 and 44, respectively The bound 
ground states of these two nuclei are thus oblate, whereas 
in the lighter isotopes the oblate solutions corresponded 
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to secondary minima. This is the origin of the sud- 
den change in two-neutron separation energies, which 
become negative in *^Mg and ^"^Mg (S'2„=— 2.237 and 
— 1.975MeV, respectively. In the case of *^Mg, however, 
two-neutron emission should be hindered by the fact that 
the parent and daughter nuclei have dramatically differ- 
ent shapes, and, by this token, ^^Mg may still have a 
substantial half-life even though it is beyond the two- 
neutron drip line. 

The bottom panel of Fig. ^ shows the neutron and 
proton rms radii, r„ and r^. At the proton drip- line, the 
neutron rms radii are smaller than the proton rms radii, 
and then they increase with increasing neutron number. 
Around ^''■^^Mg, r„ becomes almost equal to r^, and for 
nuclei close to the neutron drip-line, r„ takes significantly 
larger values than rp. The increase of r„ is fairly linear, 
similarly as in Refs. j|j5^,^, and gives no hint of an 
existence of unusually larger neutron distributions at the 
neutron drip line (see also the discussion in Refs. |^ , ^ ). 



D. Neutron-drip-line calculations 

Having at our disposal a viable method for performing 
deformed HFB calculations up to the drip lines, we have 
performed a systematic study of the equilibrium prop- 
erties of the neutron-rich nuclei in all even-Z isotopic 
chains with proton numbers from Z—2-18. In this way, 
we have explored the neighborhood of the neutron-drip 
line for all neutron numbers from A'^=6-40. 

We first performed spherical HFB-f SLy4 calculations 
in coordinate space, using the methods and the code de- 
veloped in Ref. Q. We used volume delta pairing, with 
a coupling constant Vb=— 218.5 MeVfm'^, adjusted as in 
Ref. . This value is very c lose t o the one used in our 
deformed THO code (see Sec. [V A ), suggesting that the 
effective pairing phase spaces used in the two approaches 
are very similar to one another. 

From the spherical calculations, we obtained that the 
heaviest even isotopes, for which the Fermi energies are 
negative are: ^He, ^^B, 22c, 289, ^o^e, '^"Mg, ''^^Si, ^°S, 
^*Ar. We used these spherical results as a starting point 
for our deformed calculations. 

Next, within the deformed THO formalism, we found 
that the heaviest isotopes with negative Fermi energies 
are: ^He, ^^B, 22c, 280^ se^g, ^^Mg, ^^Si, ^^S, ^SAr. 
The results obtained for these nuclei are summarized in 
Fig. 1^. By comparing the deformed results to the spher- 
ical results, we see that the position of the last bound 
nucleus is influenced by deformation only in ^^Ne and 
^^S. Volume pairing correlations are very weak in these 
nuclei; indeed, in all but the one case of ^^Ne, neutron 
pairing vanishes in the last bound nucleus of an isotopic 
chain. This suggests the necessity of using a surface pair- 
ing force here. Such a conclusion is supported by the fact 
that HB calculations |Q , carried out with a Gogny pair- 
ing force, give sizable neutron pairing correlations in this 



region (Note that surface pairing and Gogny pairing pro- 
duce quite similar distributions ||7|| over the single-particle 
states.) 

Since neutron pairing vanishes in ^^B, our result is 
identical to that of Ref. , namely that the SLy4 force 
does not produce ^"^B as bound, in disagreement with 
experiment pO| . Similarly, neither pairing nor deforma- 
tion effects are present in the calculated ^^O nucleus, 
and hence this nucleus remains bound (see discussion in 
Sec. [V B| ). On the other hand, the SLy4 force correctly 
describes ^He Q and ^^C as the last bound nuclei of 
their respective isotope chains [ |6l| ]. 

A remarkable result obtained in our calculations is that 
the last bound nuclei for all chains of isotopes heavier 
than oxygen have oblate ground-state shapes. In all of 
them, the mechanism for this effect is id entica l to that 
discussed for the Mg isotopes (see Sec. IV Cj) , namely 
that the neutron Fermi energy A„, as a function of the 
neutron number N, becomes positive for smaller values of 
N in the prolate ground states than it does in the oblate 
secondary minima. Therefore, in the heaviest bound iso- 
topes, the prolate states are unbound, whereas the oblate 
states continue to be bound and become the ground-state 
configurations. 



V. CONCLUSIONS 

In this paper, we have applied a local-scaling point 
transformation to the deformed three-dimensional Carte- 
sian harmonic oscillator wave functions so as to allow 
for a modification of their unphysical asymptotic prop- 
erties. In this way, we have obtained single-particle 
bases that remain infinite, discrete, and complete, but for 
which the wave functions have the asymptotic properties 
that are required by the canonical bases of Hartree-Fock- 
Bogoliubov theory. These bases preserve all the simplic- 
ity of the original harmonic-oscillator wave functions, and 
at the same time are amenable to very efficient numerical 
methods, such as Gauss-integration quadratures. They 
also allow for very simple calculations of local densities, 
which are at the core of self-consistent methods based on 
a Skyrme effective interaction. 

The axial transformed harmonic oscillator basis has 
been implemented to achieve a fast and reliable method 
for solving the HFB equations with the correct asymp- 
totic conditions. We have discussed several practical as- 
pects of the implementation, like the treatment of pairing 
correlations, and tested the convergence and accuracy. 

The formalism was developed for a general deformed 
transformed harmonic oscillator basis. Practical applica- 
tion within the configurational HFB formalism suggested 
a simplification to a purely spherical basis, however, as 
had been used in earlier calculations. We have neverthe- 
less presented the general formalism because of its possi- 
ble use in other applications. 

As a first application of this new methodology, we have 
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carried out HFB calculations using the SLy4 Skyrme 
force and a density-independent (volume) pairing force. 
The calculations were performed for the chain of even-Z 
Mg isotopes and for the light even-Z nuclei located near 
the two-neutron drip line. We have presented results for 
binding energies, quadrupolc moments, and for the pair- 
ing properties of these nuclei. 

Perhaps the most interesting outcome of our calcula- 
tions is that nuclei that are formally beyond the two- 
neutron drip line, i.e., those with negative two-neutron 
separation energies, may have tangible half lives, pro- 
vided (i) that they have localized ground states (neg- 
ative Fermi energies), and (ii) their ground-state con- 
figurations are significantly different than those of the 
(daughter) nuclei with two less neutrons. According to 
our calculations, precisely such a situation occurs in the 
chains of isotopes with Z=10, 12, 14, 16 and 18. In these 
chains, the prolate configuration becomes unbound be- 
fore (i.e., for a smaller neutron number) than the oblate 
configuration. That change in the ground state struc- 
ture leads to negative two-neutron separation energies 
and thus to the exotic conditions given above. 
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FIG. 1. Total energies E, and proton and neutron rms radii, Vp and Vn, obtained in the HFB+SLy4 calculations for ^^0 by 
using the HO and THO bases, as functions of the number of HO shells iVgh. The exact results refer to those obtained from 
spherical coordinate-space calculations. 
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FIG. 2. Neutron densities obtained in the HFB+SLy4 calculations for by using the HO (dashed line) and THO (solid line) 
bases. Neutron and proton densities denoted as "exact" (dots) have been obtained from spherical coordinate-space calculations 
in a box of i?box =20 fm. 
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FIG. 3. Total energies E, proton and neutron rrris radii, and r„, and deformations j3 obtained in the HFB+SLy4 calcula- 
tions for ""^Mg by using the HO and THO bases, as functions of the number of HO shells Ngh. The horizontal lines denote the 
THO results obtained at iVsh=20. 
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FIG. 4. Neutron densities obtained in the HFB+SLy4 calculations for the deformed ground state of '*'*Mg by using the HO 
(circles) and THO (squares) bases. Each point corresponds to one Gauss-integration node in the z-p plane, and the results are 
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FIG. 5. Neutron Fermi energies A„, energies per particle E/A, deformations /3, quadrupole moments Q, pairing gaps A, and 
pairing energies Spair calculated for the Mg isotopes within the HFB+SLy4 method in the THO basis {Nsh='20), as functions of 
the mass number A. Apart from the upper panel, circles, squares, and diamonds pertain to proton, neutron, and total results, 
respectively. Closed symbols connected with lines denote values for the absolute minima in the deformation-energy curve (axial 
shapes are assumed), while open symbols pertain to secondary minima. 
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FIG. 6. Ujjpcr panel: two-neutron separation energies S2n (open symbols) compared to — 2A„ (closed symbols), and lower 
panel: proton and neutron rms radii. Calculations for the Mg isotopes were performed within the HFB-|-SLy4 method in the 
THO basis for iV3h=20. 
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FIG. 7. Neutron Fermi energies \n, energies per particle E/A, pairing gaps A, pairing energies E, deformations /?, quadrupole 
moments Q, and rms radii r calculated for the neutron-drip- line nuclei (indicated in the lower panel) within the HFB+SLy4 
method in the THO basis, as functions of the mass number A. Circles, squares, and diamonds pertain to proton, neutron, and 
total results, respectively. 
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